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We employed a fully optimized Shadow Wave Function (SWF) in combination with Variational 
Monte Carlo techniques to investigate the properties of HD molecules and molecular ortho-deuterium 
(0-D2) in bulk solid para- hydrogen (p-Eb). Calculations were performed for different concentrations 
of impurities ranging from about 1% to 25% at the equilibrium density for the para-hydrogen crystal. 
By computing the excess energy both for clustered and isolated impurities we tried to determine a 
limit for the solubility of HD and 0-D2 in p-Eb. 
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Isotopic impurities in quantum crystals are often used 
as a marker for testing the occurrence of quantum me- 
chanical particle-particle exchange by means of the NMR 
experiments [l|, 0]. In particular, in the case of J = 
impurities, like HD or (0-D2) molecules, embedded in a 
J = para- hydrogen (P-H2) matrix with small concen- 
trations of ortho- hydrogen (less than 3%), the observed 
longer time relaxation can only be attributed to particle- 
particle tunneling, because the so-called "resonant con- 
version" is not allowed (see Ref. and its inner refer- 
ences). Concentrations of HD molecules in these exper- 
iments are about l-2%[l|. It is safe to assume that the 
interaction among different hydrogen isotopes are very 
similar. Only the mass difference should distinguish the 
behavior of P-H2 and HD. This is the main reason for re- 
lating the motion of HD impurities to the motion of the 
molecules in the matrix. However, in a quantum mechan- 
ical solid, the mass difference among the constituents can 
act in quite subtle way. In fact, heavier particles tend to 
be more localized, with a net gain in energy, while their 
nearest neighbors might find convenient to lower their ki- 
netic energy by relaxing towards the impurity, in a very 
delicate balance between the gain in kinetic energy and 
the loss in potential energy. This is the case of HD and o- 
D2 impurities in H2. The presence of impurities of lower 
mass (like 3 He in 4 He) should instead have a less dra- 
matic effect, due to steric hindrance which prevents the 
impurity from relaxing too much against the rest of the 
crystal. Change of force constants and lattice distortion 
in the proximity of a isotopic impurity were studied in 
relation to their influence on the thermal conductivity 
of solid p-H 2 at low temperatures^. More recently, Ra- 
man spectroscopy experiments [3] clearly showed an effect 
on both rotational and vibrational spectra which can be 
attributed to local relaxation of the lattice around the 
impurity. The balance among kinetic and potential en- 
ergy might also lead to clustering phenomena when the 
concentration of impurities exceeds a certain limit. Fi- 
nally, HD impurities might play a role in the determining 
the results on measurements of non classical torsional in- 



ertia recently performed to investigate the existence of a 
supersolid phase in P-H2 1| . 

In this letter we report the results of Variational Monte 
Carlo (VMC) simulations performed with a modified 
Shadow Wave Function[6j, |7| (SWF) to study the prop- 
erties of HD and 0-D2 molecules in a pure P-H2 matrix. 
The main property of SWF is the capability of describ- 
ing both the crystalline and the disordered phase of a 
quantum system with the same functional form, with is 
always translationally invariant. Unlike in the standard 
Jastrow-Nosanow trial wave function, using SWF no a 
priori equilibrium are imposed to describe the crystalline 
phase. This property allows SWF to provide a realistic 
description of a variety of inhomogeneous system, like 
liquid-solid phase coexistence in 4 He[8|, as well as finite 
system as, for instance, 4 He clusters with and without 
impurities 0. Recently a SWF has been employed also 
for P-H2, in order to describe the homogeneous liquid 
and solid phase [lfj], as well as to investigate energetic 
and structural properties of defective crystals 11 1. 

We model the system of Njj molecules of P-H2 and 
Nj isotopic impurities as a set of N = Nh + Nj point 
particles described by the following Hamiltonian: 
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where v(r) is the well-known Silvera-Goldman (SG) 
potential[12j|. This model interaction includes a term 
oc 1/r 9 which effectively accounts for triple dipole inter- 
actions in the system. This approximation is supported 
by the fact that both the P-H2 molecule and the isotopic 
impurities considered in this paper, HD and o-D 2 are in 
the rotational ground state ( J = 0), and therefore the av- 
erage interaction must be spherically symmetric. Higher 
terms in a multipolar expansion, due to the non-spherical 
nature of the molecules, may be safely considered as sec- 
ond order. The limit of validity of the SG potential and 
the influence of explicit many-body terms in the poten- 
tial on the equation of state have been recently tested 
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both for solid p-EbQjl and 0-D2 14 1. In the latter case, 
the SG potential has been found to give a rather good 
description of the equation of state, while in the case of 
solid p-H 2 the same interaction leads to a slightly more 
inaccurate result on the energies due to stronger many- 
body effects. However, the explicit angular dependence 
of the potential seems not to have sizeable effects on the 
structural properties, such as the pair distribution func- 
tion. The many-body quantum mechanical problem is 
solved at temperature T=0 by means of a trial solution 
for the ground-state of the Hamiltonian in Eq. [T] in the 
form of a modified Shadow Wave Function: 
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where R = {ri...r/v} are the coordinates of the 
molecules and S = {si . . . Sjv} are 37V auxiliary degrees of 
freedom named shadows. Both functions (f> r and 4> s have 
been expressed as Jastrow products of pair functions: 
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while the kernel K (R, S) is written as a product of Gaus- 
sians connecting the real and auxiliary degrees of free- 
dom: 
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The pair pseudopotentials between molecules have been 
expanded in term of suitable basis functions: 
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where both x and y indicate generically one of the two 
label " H" or " 7" . The basis functions Xn (about 40) are 
the same used in Ref. [l3j and [lj for DMC calculations 
in solid P-H2 and 0-D2. 
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where L is the side of the simulation box. The cutoff 
radius r c allows to remove divergences for r — > in or- 
der to avoid numerical instabilities in the optimization 
procedure. If such cutoff is small enough, it does not in- 
fluence the value of the energy. In our calculations it has 
been taken equal to lA. The pseudopotential between 



FIG. 1: Excess energy per impurity (in K) in the P-H2 crystal 
as function of the number of 0-D2 impurities. Filled triangles: 
close impurities, empty circles: far impurities. 



auxiliary degrees of freedom is indeed just the rescaled 
Silvera-Goldman interaction: 



*Uxx(s) — fixx'vi.&xxS) • 
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The parameters {b xy }, {a™ }, {C x }, {a xy }, {5 xy } have 
been optimized following the reweighting scheme, which 
alternate between VMC simulations and minimization of 
a combination of the variance on the expectation value 
of the Hamiltonian, and of the expectation value itself 
estimated on a subset of the sampled configurations. We 
performed simulations for a face centered cubic (fee) lat- 
tice at densities p = 0.02609A~ 3 , very close to the bulk 
equilibrium density of pure p-H 2 [13| . The simulation cell 
was set to accommodate 3x3x3 elementary cubic cells 
for fee, with a total of iV=108 lattice sites. All the in- 
termolecular interactions are truncated at the edge of a 
sphere of radius equal to L/2, where L is the length of 
the side of the cell. The contribution from the potential 
energy outside the sphere is estimated by integrating the 
potential in the {L/2, +00) interval, assuming therefore 
all the pair correlation functions of molecules and im- 
purities to be constant beyond L/2. Calculations were 
performed for different number of impurities, from 1 to 
27, corresponding to concentrations ranging from about 
1% to 25%. For each value of concentration, we repeated 
the optimization of the trial wave function both with all 
the impurities initially localized around next-neighbours 
site, as close to each other as possible, and with far im- 
purities, i. e., separated by at least a shell of P-H2. The 
excess energy per impurity ei mp is defined as follow: 



Eimp = -rr [E(N H) Ni) 
iVj 



E(N,0)}, 



(8) 



where E{Nh , N{) is the total energy of the system of Njj 
molecules of p-H 2 and Nj = N — Nh impurities. It can 
be noticed that the excess energy per impurity remains 
almost constant in the range of concentrations consid- 
ered. A rough estimate of the excess energy of a single 
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P-H2 close 


P-H2 far 


O-D2 close 


O-D2 far 


V 


-158,15(2) 


-160,67(2) 


-171.49(3) 


-164.71(3) 


Ekin 


69,76(3) 


71.15(3) 


47.82(4) 


45.57(3) 


V + Ekin 


-88,39(2) 


-89,52(3) 


-123.67(3) 


-119.14(3) 



TABLE I: Potential, kinetic and total energy (in Kelvin) per 
particle of p-Eb and 0-D2 molecules for a p-Irb crystal with 
25% of 0-D2 impurities. Labels "far" and "close" indicate sys- 
tems with diluted or clustered impurities respectively. 



4 8 12 16 20 24 28 

N i 

FIG. 2: Excess energy per impurity (in K) in the P-H2 crystal 
as function of the number of HD impurities. Filled triangles: 
close impurities, empty circles: far impurities. 

isotopic impurity in the p-H 2 crystal might be obtained 
considering as a perturbation term in the Hamiltonian 
the difference in kinetic energy AEkin due to the mass 
difference of the two isotopes. 

&imp ~ ( 1 ) A-^kin (9) 

\m D J 

The excess binding energy estimated in this way is about 
35K for 0-D2 and about 24K for HD. Both values are 
not far from the values of about 40K and 25K respec- 
tively, obtained from the VMC simulations. This fact 
suggests that the contribution of the mass difference to 
the excess energy of the impurity is dominant, in par- 
ticular for HD impurities. We point out that in SWF 
simulations both relaxation and possible quantum diffu- 
sion effects are correctly taken into account. As shown 
in Figures [T] and [2] the excess energy computed with 
clustered or diluted impurities, is essentially coincident 
within errorbars for low concentrations. For concentra- 
tions up to 10% we are therefore lead to consider the 0-D2 
as not subject to clustering in the P-H2 crystal matrix. 
For higher concentrations, however, the excess energy per 
impurity for clustered molecules becomes lower than the 
energy for the diluted ones well outside errorbars. In 
the case of HD impurities, this threshold is higher, as it 
might be expected because of the lighter mass, and the 
VMC results are compatible with a limiting concentra- 
tion for clustering of about 15%. In order to check if the 
estimate of the excess energy and of the clustering limits 
are affected either by finite-size effects or by the choice 
of the lattice fitted by the simulation box, we performed 
some simulations in a larger fee lattice using N = 256 
molecules arranged on 4x4x4 elementary cubic cells and 
in the hep lattice using N = 180 molecules filling 5x3x3 
elementary cells. In this case the simulation box is not 
cubic, but the ratio of the three sides is as close as pos- 
sible to 1. The results for concentration of 25% of 0-D2 
molecules obtained with different lattices and different 





P-H2 close 


P-H2 far 


HD close 


HD far 


V 


-158,05(2) 


-159,63(2) 


-167.03(3) 


-161.78(3) 


Ekin 


70,01(3) 


70.73(3) 


57.51(4) 


55.16(4) 


V + Ekin 


-88,04(2) 


-88,90(2) 


-109.52(3) 


-106.62(3) 



TABLE II: Potential, kinetic and total energy (in Kelvin) per 
particle of P-H2 and HD molecules for a P-H2 crystal with 25% 
of HD impurities. Labels "far" and "close" indicate systems 
with diluted or clustered impurities respectively. 



values of N show that the values of the excess energy 
are largely independent from from the total number of 
particles, indicating the absence of substantial size-finite 
effects. 

The energy differences between the system with clus- 
tered or diluted impurities, can be splitted into the 
mean kinetic and potential energy contributions for p- 
H2 molecules and impurities, as reported in Tables [J and 
ITT1 From this analysis it is possible to see how the lower 
energy of the system with clustered impurities is due to 
the drastic lowering of the expectation of the potential 
energy of the impurities themselves. As shown below, 
this result is related to the stronger localization of the 
impurities around equilibrium positions on a lattice with 
spacing reduced compared to that of P-H2 molecules. Al- 
though the average total energy of the p-H 2 molecules 
becomes higher when the impurities are clustered, the 
net effect is lower total energy of the system. Some in- 
formation about the distortions in the crystal induced 
by the impurities in the P-H2 matrix can be obtained by 
computing the pair correlation function 

<Mr) = ^££(|r?-rf|-r)^, (10) 

where a and /3 indicate the species among which dis- 
tances are computed and from the configurations of the 
system sampled in the MC simulations. When using the 
SWF it is also possible to define the same operator for 
the shadow degree of freedom 

g sa p(s) = (j2s(K-s^\-s)\, (11) 

As discussed elsewhere [H, the shadows in quantum 
crystals are more rigid degrees of freedom, which act 
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FIG. 3: Pair correlations functions g a /3(r) and g S a/3(r) for 
a = f3 =p-Hh and a — /3=o-D2 for a 25% concentration of 
clustered impurities in a fee lattice. Solid and dashed arrows 
indicate the positions of the next-neighbours peaks for p-H2 
and 0-D2 molecules, respectively. 
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FIG. 4: Pair correlations functions g a p(r) and g sa f3(r) with 
a — (3 =p-H2 and a = (5 =HD for a 25% concentration of 
clustered impurities in a fee lattice. Solid and dashed arrows 
indicate the positions of the next-neighbours peaks for P-H2 
and HD molecules, respectively. 

as heavier particles indicating a sort of "average" po- 
sition for the real quantum particles. In Figures [3] and 
[4] we report the pair correlations functions both for P-H2 
molecules and for the impurities in the case of a con- 
centration of 25% of clustered impurities embedded in a 
fee P-H2 matrix. The 0-D2 impurities tend to be more 
localized than the P-H2 molecules, as expected due the 
higher mass. Moreover, the peaks of the pair correlation 
function of clustered 0-D2 impurities are clearly shifted 
inwards, with respect to the pair correlation of surround- 
ing hydrogen molecules. The strong localization of the 
impurities around equilibrium positions closer to each 
other than ones of P-H2 give probably rise to the drastic 
lowering of the potential energy of the 0-D2 molecules 
discussed above, making the clustering of impurities en- 



ergetically favourite. 

In conclusion, we used a Shadow Wave Function in 
connection with Variational Monte Carlo calculations to 
study the properties of 0-D2 and HD molecules in a p- 
H2 crystal. In particular we computed the excess energy 
for clustered and isolated impurities at different concen- 
trations, trying to determine a limit for the solubility of 
o-D 2 and HD in p-H 2 . 

For concentrations up to 10% and 15%, for 0-D2 and 
HD molecules respectively, the excess energy computed 
for clustered impurities becomes lower than that for di- 
luted ones. 

Calculations of the pair correlation function show that 
at large concentrations clustered impurities are localized 
around equilibrium positions closer than ones of P-H2 
molecule giving rise to a drastic lowering of the potential 
energy of the impurities themselves. 

We are grateful to C.J. Umrigar and P. Nightingale 
for providing us the Levemberg-Marquardt minimization 
code used for optimizing the parameters in the varia- 
tional wave function. 

Calculations have been performed partly on the HPC 
facility of the Departement of Physics, University of 
Trento, and partly at the computing facilities at ECT* 
in Trento. 
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